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ABSTRACT 

The Black-Scholes theory of option pricing has been considered for many years as an 
important but very approximate zeroth-order description of actual market behavior. We 
generalize the functional form of the diffusion of these systems and also consider multi-factor 
models including stochastic volatility. Daily Eurodollar futures prices and implied volatilities 
are fit to determine exponents of functional behavior of diffusions using methods of global 
optimization, Adaptive Simulated Annealing (ASA), to generate tight fits across moving time 
windows of Eurodollar contracts. These short-time fitted distributions are then developed 
into long-time distributions using a robust non-Monte Carlo path-integral algorithm, 
PATHINT, to generate prices and derivatives commonly used by option traders. 
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1. INTRODUCTION 

1.1. Background 

There always is much interest in developing more sophisticated pricing models for financial instruments. 
In particular, there currently is much interest in improving option pricing models, particularly with respect 
to stochastic variables [1-4]. 

The standard Black-Scholes (BS) theory assumes a lognormal distribution of market prices, i.e., a 
diffusion linearly proportional to the market price. However, many texts include outlines of more general 
diffusions proportional to an arbitrary power of the market price [5]. 

The above aspects of stochastic volatility and of more general functional dependencies of diffusions 
are most often "swept under the rug" of a simple lognormal form. Experienced traders often use their 
own intuition to put volatility "smiles" into the BS theoretical constant coefficient in the BS lognormal 
distribution to compensate for these aspects. 

It is generally acknowledged that since the market crash of 1987, markets have been increasingly 
difficult to describe using the BS model, and so better modelling and computational techniques should be 
used traders [6], although in practice simple BS models are the rule rather than the exception simply 
because they are easy to use [7]. To a large extent, previous modelling that has included stochastic 
volatility and multiple factors has been driven more by the desire to either delve into mathematics 
tangential to these issues, or to deal only with models that can accommodate closed-form algebraic 
expressions. We do not see much of the philosophy in the literature that has long driven the natural 
sciences: to respect first raw data, secondly models of raw data, and finally the use of numerical 
techniques that do not excessively distort models for the sake of ease of analysis and speed of 
computation. Indeed, very often the reverse set of priorities is seen in mathematical finance. 

1.2. Our Approach 

We have addressed the above issues in detail within the framework of a previously developed 
statistical mechanics of financial markets (SMFM) [8-13]. 

Our approach requires three sensible parts. Part one is the formulation of the model, which to some 
extent also involves specification of the specific market(s) data to be addressed. Part two is the fitting of 
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the model to specific market data. Part three is the use of the resulting model to calculate option prices 
and their Greeks (partial derivatives of the prices with respect to their independent variables), which are 
used as risk parameters by traders. Each part requires some specific numerical tuning to the market under 
consideration. 

The first part was to develop the algebraic model to replace/generalize BS, including the possibility 
of also addressing how to handle data regions not previously observed in trading. This is not absurd; 
current BS models perform integrals that must include a much influence from fat tails that include data 
regions never seen or likely to be seen in real-world markets. There are some issues as to whether we 
should take seriously the notion that the market is strongly driven by some element of a "self-fulfilling 
prophesy" by the BS model [14], but in any case our models have parameters to handle a wide range of 
possible cases that might arise. 

We have developed two parallel tracks starting with part one, a one-factor and a two-factor model. 
The two-factor model includes stochastic volatility. At first we sensed the need to develop this two-factor 
model, and we now see that this is at the least an important benchmark against which to judge the worth 
of the one-factor model. 

The second part was to fit the actual raw data so we can come up with real distributions. Some tests 
illustrated that standard quasi-linear fitting routines, could not get the proper fits, and so we used a more 
powerful global optimization, Adaptive Simulated Annealing (ASA) [15]. Tuning and selection of the 
time periods to perform the fits to the data were not trivial aspects of this research. Practical decisions 
had to be made on the time span of data to be fit and how to aggregate the fits to get sensible "fair values" 
for reasonable standard deviations of the exponents in the diffusions. 

The third part was to develop Greeks and risk parameters from these distributions without making 
premature approximations just to ease the analysis. Perhaps someday, simple approximations and 
intuitions similar to what traders now use for BS models will be available for these models, but we do not 
think the best approach is to start out with such approximations until we first see proper calculations, 
especially in this uncharted territory. When it seemed that Cox-Ross-Rubenstein (CRR) standard tree 
codes (discretized approximations to partial differential equations) [16] were not stable for general 
exponents, i.e., for other than the lognormal case, we turned to a PATHINT code developed a decade ago 
for some hard nonlinear multifactor problems [17], e.g., combat analyses [18], neuroscience [19,20], and 
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potentially chaotic systems [21,22]. In 1990 and 1991 papers on financial applications, it was mentioned 
how these techniques could be used for stochastic interest rates and bonds [9,10]. The modifications 
required here for one-factor European and then American cases went surprisingly smoothly; we still had 
to tune the meshes, etc. The two-factor model presented a technical problem to the algorithm, which we 
have reasonably handled using a combination of selection of the model in part one and a reasonable 
approach to developing the meshes. 

1.3. Outline of Paper 

Section 1 is this introduction. Section 2 describes the nature of Eurodollar (ED) futures data and 
the evidence for stochastic volatility. Section 3 outlines the algebra of modelling options, including the 
standard BS theory and our generalizations. Section 4 outlines the three equivalent mathematical 
representations used by SMFM; this is required to understand the development of the short-time 
distribution that defines the cost function we derive for global optimization, as well as the numerical 
methods we have developed to calculate the long-time evolution of these short-time distributions. Section 
5 outlines ASA and explains its use to fit short-time probability distributions defined by our models to the 
Eurodollar data; we offer the fitted exponent in the diffusion as a new important technical indicator of 
market behavior. Section 6 outlines PATHINT and explains its use to develop long-time probability 
distributions from the fitted short-time probability distributions, for both the one-factor and two-factor 
tracks. Section 7 describes how we use these long-time probability distributions to calculate European 
and American option prices and Greeks; here we give numerical tests of our approach to BS CRR 
algorithms. Section 8 is our conclusion. 

2. DATA 

2.1. Eurodollars 

Eurodollars are fixed-rate time deposits held primarily by overseas banks, but denominated in US 
dollars. They are not subject to US banking regulations and therefore tend to have a tighter bid-ask 
spread than deposits held in the United States [23]. 
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2.2. Futures 

The three-month Eurodollar futures contract is one of the most actively traded futures markets in 
the world. The contract is quoted as an index where the yield is equal to the Eurodollar price subtracted 
from 100. This yield is equal to the fixed rate of interest paid by Eurodollar time deposits upon maturity 
and is expressed as an annualized interest rate based on a 360-day year. The Eurodollar futures are cash 
settled based on the 90-day London Interbank Offer Rate (LIBOR). A "notional" principal amount of $1 
million, is used to determine the change in the total interest payable on a hypothetical underlying time 
deposit, but is never actually paid or received [23]. 

Currently a total of 40 quarterly Eurodollar futures contracts (or ten years worth) are listed, with 
expirations annually in March, June, September and December. 

2.3. Options on Futures 

The options traded on the Eurodollar futures include not only 1 8 months of options expiring at the 
same time as the underlying future, but also various short dated options which themselves expire up to 
one year prior to the expiration of the underlying futures contract. 

2.4. Front/Back Month Contracts 

For purposes of risk minimization, as discussed in a previous paper [4], traders put on spreads 
across a variety of option contracts. One common example is to trade the spread on contracts expiring 
one year apart, where the future closer to expiration is referred to as the front month contract, and the 
future expiring one year later is called the back month. The availability of short dated or "mid-curve" 
options which are based on an underlying back month futures contract, but expire at the same time as the 
front month, allow one to trade the volatility ratios of the front and back month futures contracts without 
having to take the time differences in option expirations into consideration. We studied the volatilities of 
these types of front and back month contracts. Here, we give analyses with respect only to quarterly data 
longer than six months from expiration. 
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2.5. Stochastic Volatility 

Below we develop two-factor models to address stochastic volatility. In a previous paper, we have 
performed empirical studies of Eurodollar futures to support the necessity of dealing with these issues [4]. 

3. MODELS 

3.1. Random walk model 

The use of Brownian motion as a model for financial systems is generally attributed to 
Bachelier [24], though he incorrectly intuited that the noise scaled linearly instead of as the square root 
relative to the random log-price variable. Einstein is generally credited with using the correct 
mathematical description in a larger physical context of statistical systems. However, several studies 
imply that changing prices of many markets do not follow a random walk, that they may have long-term 
dependences in price correlations, and that they may not be efficient in quickly arbitraging new 
information [25-27]. A random walk for returns, rate of change of prices over prices, is described by a 
Langevin equation with simple additive noise ?], typically representing the continual random influx of 
information into the market. 

f = -n + rm , 

f = dTldt , 

< r](t) > v =0,< rj(t), rj(f) > n = 8{t - t') , (1) 

where y\ and y 2 are constants, and T is the logarithm of (scaled) price. Price, although the most dramatic 
observable, may not be the only appropriate dependent variable or order parameter for the system of 
markets [28]. This possibility has also been called the "semistrong form of the efficient market 
hypothesis" [25]. 

The generalization of this approach to include multivariate nonlinear nonequilibrium markets led to 
a model of statistical mechanics of financial markets (SMFM) [8]. 
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3.2. Black-Scholes (BS) Theory 

The standard partial-differential equation used to formulate most variants of Black-Scholes (BS) 
models describing the market value of an option, V, is 

dV 1 , , d 2 V dV 

where S is the asset price, and a is the standard deviation, or volatility of S, and r is the short-term 
interest rate. The solution depends on boundary conditions, subject to a number of interpretations, some 
requiring minor transformations of the basic BS equation or its solution. For example, the basic equation 
can apply to a number of one-dimensional models of interpretations of prices given to V, e.g., puts or 
calls, and to S, e.g., stocks or futures, dividends, etc. 

For instance, if V is set to C, a call on an European option with exercise price X with maturity at T, 
the solution is 

C(S, t) = SN{d x )- Xe-^^Nidi) , 
ln(S/X) + (r+ l ~ o 2 )(T - t) 

dl= <7(J-0 1/2 ' 



ln(5/X) + (r - - cr 2 )(r - t) 
di = ^Tn • (3) 

In practice, the volatility a is the least known parameter in this equation, and its estimation is 
generally the most important part of pricing options. Usually the volatility is given in a yearly basis, 
baselined to some standard, e.g., 252 trading days per year, or 360 or 365 calendar days. Therefore, all 
values of volatility given in the graphs in this paper, based on daily data, would be annualized by 



multiplying the standard deviations of the yields by V252 = 15. 87. We have used this factor to present our 
implied volatilities as daily movements. 



3.3. Some Key Issues in Derivation of BS 

The basic BS model considers a portfolio in terms of delta (A), 

n = v-as , 



(4) 
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in a market with Gaussian-Markovian ("white") noise X and drift ju, 
dS 

— = adX + udt , 
S 

where V(S, t) inherits a random process from S, 



(5) 



dV 

dV = aS — dX + 
dS 



'dV 1 2 ,d 2 V dV^ 
juS — + - a 2 S 2 — + 



dS 2 dS 2 dt 



dt . 



(6) 



This yields 



J 



dV 1 2 c 2 d2y dV 



dS 2 

The expected risk-neutral return of II is 
dU = rUdt = r(V-AS)dt . 



juS — + - a^S^ — y + ^— - juAS dt 



dS 2 dt 



(7) 



(8) 



Options V on futures F can be derived, e.g., using simple transformations to take cost of carry into 
consideration, such as 



F = Se r(T - f) , 
and setting 

dU = rV dt . 
The corresponding BS equation for futures F is 
d 2 V 



(9) 



(10) 



Ik 



+ -a 2 F 2 



dS 2 



-rV = 



(11) 



At least two advantages are present if A is chosen such that 
dV 



A = 



dS 



(12) 



Then, the portfolio can be instantaneously "risk-neutral," in terms of zeroing the coefficient of X, as well 
as independent of the direction of market, in terms of zeroing the coefficient of ju. For the above example 
of V = C, 



A=N(d l ) 



(13) 
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Other trading strategies based on this simple model use similar constructs as risk parameters, e.g., 
gamma (T), theta (0), vega (T), rho (p) [5], 

a 2 n 



r = 



as 2 ' 



an 



T = 



an 

a^ ' 

an 

dr 



The BS equation, Eq. (2), may be written as 



(14) 



® + rSA+^(aSyr=rf 



(15) 



3.4. S x Models 

Our two-factor model includes stochastic volatility a of the underlying 5, 

dS = pdt + a F(S, S , S^, x, y) dzs , 

da = v dt + e dz a , 

< dn >= , i = {S, a} , 



< dn(t) dzj(t') > = dt S(t - 1') ,i = j , 



< dzi(t) dzj{f) > = p dt 8{t - f) ,i*j, 



F(S, S , S^, x, y) = 



s, 

cx cl—x 



s<s 
s Q <s<s c 



(16) 



syst x s x - y , s>s c 



where Sq and 5^ are selected to lie outside the data region used to fit the other parameters, e.g., S = 1 
and = 20 for fits to Eurodollar futures which historically have a very tight range relative to other 
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markets. We have used the Black-Scholes form F = S inside 5 < S to obtain the usual benefits, e.g., no 
negative prices as the distribution is naturally excluded from 5 < and preservation of put-call parity. 
Put-call parity for European options is derived quite independent of any mathematical model of 
options [5]. In its simplest form, it is given by 



c + Xe- r(T - t] = p + S , 



(17) 



where c (p) is the fair price of a call (put), X is the strike price, r is the risk-free interest rate, t is the 
present time, T is the time of expiration, and S is the underlying market. We have taken y = 0, a normal 
distribution, to reflect total ignorance of markets outside the range of 5 > 5^. The one-factor model just 
assumes a constant a. It is often noted that BS models incorrectly include untenable contributions from 
large 5 regions because of their fat tails [29]. (If we wished to handle negative interest rates, ED prices > 
100, we would move shift the 5 = axis to some 5 < value.) 

We found that the abrupt, albeit continuous, changes across 5 especially for x < did not cause 
any similar effects in the distributions evolved using these diffusions, as reported below. 

The formula for pricing an option P, derived in a Black-Scholes generalized framework after 
factoring out interest-rate discounting, is equivalent to using the form 

dS = ju S dt + a F(S, 5 , S M , x, y) dis , 



da = v dt + e dz^ 



(18) 



We experimented with some alternative functional forms, primarily to apply some smooth cutoffs 
across the above three regions of 5. For example, we used F', a function F designed to revert to the 
lognormal Black-Scholes model in several limits, 

F'(S, So, Sw x) = 5 C + (1 - C ) (S x + 5 (1 - C^)) , 



Cq = exp 



Coo = exp 



fs_ ii — jci y-* i+r 

v 5 1 + ll-jcl, 



(— T 

Woo/ 



lim F'(S, So, Soo, x) = Sq = constant , 

S—>oo, X*\ 
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lim F'(S, S , S.., x) = lim F'(S, S Q , S„, x) = S . (19) 

However, our fits were most sensitive to the data when we permitted the central region to be pure S x using 
F above. 

3.4.1. Various F(S, x) Diffusions 

Fig. 1 gives examples of F(S,So,S oa ,x,y)dzs for x in {-1, 0, 1, 2}. The other parameters are 
5 = 5, S = 0. 5, 5^ = 20, y = 0. 



Fig. 1. 



4. STATISTICAL MECHANICS OF FINANCIAL MARKETS (SMFM) 

4.1. Statistical Mechanics of Large Systems 

Aggregation problems in nonlinear nonequilibrium systems typically are "solved" (accommodated) 
by having new entities/languages developed at these disparate scales in order to efficiently pass 
information back and forth. This is quite different from the nature of quasi-equilibrium quasi-linear 
systems, where thermodynamic or cybernetic approaches are possible. These approaches typically fail for 
nonequilibrium nonlinear systems. 

Many systems are aptly modeled in terms of multivariate differential rate-equations, known as 
Langevin equations, 

M° = f + gj?] j , (G = 1, • • • , A) , ( j = 1, • • • , AO , 
M G = dM G ldt , 

< 7] j (t) >ri = , < rj j {t), TjJ'it') >„= S jj 'S(t - O , (20) 

where f G and gj are generally nonlinear functions of mesoscopic order parameters M G , j is a 
microscopic index indicating the source of fluctuations, and N > A. The Einstein convention of summing 
over repeated indices is used. Vertical bars on an index, e.g., Ijl, imply no sum is to be taken on repeated 
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indices. 

Via a somewhat lengthy, albeit instructive calculation, outlined in several other papers [8,10,30], 
involving an intermediate derivation of a corresponding Fokker-Planck or Schrodinger-type equation for 
the conditional probability distribution P[M(t)\M(t )], the Langevin rate Eq. (20) is developed into the 
more useful probability distribution for M G at long-time macroscopic time event t = (u + \)9 + 1 , in 
terms of a Stratonovich path-integral over mesoscopic Gaussian conditional probabilities [31-35]. Here, 
macroscopic variables are defined as the long-time limit of the evolving mesoscopic system. 

The corresponding Schrodinger-type equation is [33,34] 
dPIdt = l - (g GG >), Gcr - (g G P\ G + V , 

g GG ' = k T S*g G g? , 

[• • -],g = d[- • -VdM G . (21) 

This is properly referred to as a Fokker-Planck equation when V = 0. Note that although the partial 
differential Eq. (21) contains information regarding M G as in the stochastic differential Eq. (20), all 
references to j have been properly averaged over. I.e., g G in Eq. (20) is an entity with parameters in both 
microscopic and mesoscopic spaces, but M is a purely mesoscopic variable, and this is more clearly 
reflected in Eq. (21). 

The path integral representation is given in terms of the "Feynman" Lagrangian L. 
P[M t \M t0 ]dM(t) = J • • • J DM exp(-S)S[M(t ) = M ]S[M(t) = M t ] , 

t 

S = kj l min J dt'L , 

u+l 

DM = lim Hg m ll(2xey m dM G , 

p=\ G 

L(M G , M G , t) = - (M G - h G )g GG '(M G ' - h G ') + - h G . G + R/6-V , 
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h G = g G - l - g - m ( g m g GG \ G >, 
gGG' = ig GG Y , 

g = det(g G G') , 

h ;G = h ,G +l GF h =8 ig h ),G > 

T F JK = g LF [JK, L] = g LF {g JLX + g KLJ - g jk,l) , 
R = g JL RjL = g JL g JK R F jKL , 

RfJKL = - (gFKJL ~ gjK,FL ~ gFLJK + gjL,FK> + gMN^FK^JL ~ ^FL^Jk) ■ (22) 

Mesoscopic variables have been defined as M G in the Langevin and Fokker-Planck representations, in 
terms of their development from the microscopic system labeled by j. The Riemannian curvature term R 
arises from nonlinear £gg'> which is a bona fide metric of this space [33]. Even if a stationary solution, 

■ G ' G 

i.e., M = 0, is ultimately sought, a necessarily prior stochastic treatment of M terms gives rise to these 
Riemannian "corrections." Even for a constant metric, the term h ;G contributes to L for a nonlinear 

mean h G . V may include terms such as X^t'g^ G > where the Lagrange multipliers J T ' G are constraints 

r 

on M G , which are advantageously modeled as extrinsic sources in this representation; they too may be 
time-dependent. 

For our purposes, the above Feynman Lagrangian defines a kernel of the short-time conditional 
probability distribution, in the curved space defined by the metric, in the limit of continuous time, whose 
iteration yields the solution of the previous partial differential equation Eq. (21). This differs from the 
Lagrangian which satisfies the requirement that the action is stationary to the first order in dt — the 
WKBJ approximation, but which does not include the first-order correction to the WKBJ approximation 
as does the Feynman Lagrangian. This latter Lagrangian differs from the Feynman Lagrangian, 
essentially by replacing R/6 above by R/12 [36]. In this sense, the WKBJ Lagrangian is more useful for 
some theoretical discussions [37]. However, the use of the Feynman Lagrangian coincides with the 
numerical method we present here using the PATHINT code. 
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Using the variational principle, J TG may also be used to constrain M G to regions where they are 
empirically bound. More complicated constraints may be affixed to L using methods of optimal control 
theory [38]. With respect to a steady state P, when it exists, the information gain in state P is defined by 

T[P] = J... j DM' Pln(P/P) , 

DM' = DM/dM u+l . (23) 

In the economics literature, there appears to be sentiment to define Eq. (20) by the Ito, rather than 
the Stratonovich prescription. It is true that Ito integrals have Martingale properties not possessed by 
Stratonovich integrals [39] which leads to risk-neural theorems for markets [40,41], but the nature of the 
proper mathematics should eventually be determined by proper aggregation of relatively microscopic 
models of markets. It should be noted that virtually all investigations of other physical systems, which are 
also continuous time models of discrete processes, conclude that the Stratonovich interpretation coincides 
with reality, when multiplicative noise with zero correlation time, modeled in terms of white noise tj j , is 
properly considered as the limit of real noise with finite correlation time [42]. The path integral 
succinctly demonstrates the difference between the two: The Ito prescription corresponds to the prepoint 
discretization of L, wherein 6M{t) — > M p+l - M p and M{t) — > M p . The Stratonovich prescription 
corresponds to the midpoint discretization of L, wherein 6M{t) — > M p+i - M p and 

M(t) — > ^ {M p+ i + M p ). In terms of the functions appearing in the Fokker-Planck Eq. (21), the Ito 

prescription of the prepoint discretized Lagrangian, L 7 , is relatively simple, albeit deceptively so because 
of its nonstandard calculus. 

L 7 (M G , M G , t)= 1 - (M G - g G )g GG iM G ' - g G ) - V . (24) 

In the absence of a nonphenomenological microscopic theory, the difference between a Ito prescription 
and a Stratonovich prescription is simply a transformed drift [36]. 

There are several other advantages to Eq. (22) over Eq. (20). Extrema and most probable states of 
are simply derived by a variational principle, similar to conditions sought in previous 
studies [43]. In the Stratonovich prescription, necessary, albeit not sufficient, conditions are given by 

<>gL = L G - L G . t = , 
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L 6:t = L 6G >M G ' + L 6& M G ' . (25) 

' G ~ ~ G ~ G 

For stationary states, M =0, and dLldM = defines « M », where the bars identify stationary 
variables; in this case, the macroscopic variables are equal to their mesoscopic counterparts. [Note that L 
is not the stationary solution of the system, e.g., to Eq. (21) with dP/dt = 0. However, in some cases [44], 
L is a definite aid to finding such stationary states.] Many times only properties of stationary states are 

" G 

examined, but here a temporal dependence is included. E.g., the M terms in L permit steady states and 
their fluctuations to be investigated in a nonequilibrium context. Note that Eq. (25) must be derived from 
the path integral, Eq. (22), which is at least one reason to justify its development. 



4.2. Correlations 

Correlations between variables are modeled explicitly in the Lagrangian as a parameter usually 
designated p (not to be confused with the Rho Greek calculated for options). This section uses a simple 
two-factor model to develop the correspondence between the correlation p in the Lagrangian and that 
among the commonly written Weiner distributions dz. 

Consider coupled stochastic differential equations 
dr = f r (r,l)dt + g r (r,l)a l dzi , 

dl = f\r, l)dt + g\r, l)a 2 dz 2 , 

<dzi >= 0,i = {1,2}, 

< dzi(t)dzj(t') >= dtS(t - 1') ,i = j, 

< dzi(t)dzj(t') >= pdtS(t - O , i 5fc j , 

«-«-$:\-,\ 

where < . > denotes expectations. 

These can be rewritten as Langevin equations (in the Ito prepoint discretization) 

drldt = f + g r o-i(y + n l + sgnp y~n 2 ) , 
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dlldt = g l + g l a 2 (sgap y n l + y + n 2 ) , 



Hi = {dt) m Pi , 

where p x and p 2 are independent [0,1] Gaussian distributions. 

The equivalent short-time probability distribution, P, for the above set of equations is 

P = g m {2ndty m exp(-Ldt) , 



(27) 



L = -F"g_F , 



F = 



(drldt-f) 
{ dlldt - f) 



g = det(g) , 



k = l-p 2 



(28) 



g, the metric in {r, Z}-space, is the inverse of the covariance matrix, 



g = 



Pg r g l O\Oi (g l 0i) 2 



(29) 



The above also corrects previous papers which inadvertently dropped the sgn factors in the 
above [9,10,17]. 



5. ADAPTIVE SIMULATED ANNEALING (ASA) FITS 



5.1. ASA Outline 

The algorithm developed which is now called Adaptive Simulated Annealing (ASA) [45] fits short- 
time probability distributions to observed data, using a maximum likelihood technique on the Lagrangian. 
This algorithm has been developed to fit observed data to a theoretical cost function over a D-dimensional 
parameter space [45], adapting for varying sensitivities of parameters during the fit. The ASA code can 



High-resolution path-integral . 



- 17- 



Lester Ingber 



be obtained at no charge, via WWW from http://www.ingber.com/ or via FTP from ftp.ingber.com [15]. 

5.1.1. General description 

Simulated annealing (SA) was developed in 1983 to deal with highly nonlinear problems [46], as an 
extension of a Monte-Carlo importance-sampling technique developed in 1953 for chemical physics 
problems. It helps to visualize the problems presented by such complex systems as a geographical terrain. 
For example, consider a mountain range, with two "parameters," e.g., along the North-South and 
East- West directions. We wish to find the lowest valley in this terrain. SA approaches this problem 
similar to using a bouncing ball that can bounce over mountains from valley to valley. We start at a high 
"temperature," where the temperature is an SA parameter that mimics the effect of a fast moving particle 
in a hot object like a hot molten metal, thereby permitting the ball to make very high bounces and being 
able to bounce over any mountain to access any valley, given enough bounces. As the temperature is 
made relatively colder, the ball cannot bounce so high, and it also can settle to become trapped in 
relatively smaller ranges of valleys. 

We imagine that our mountain range is aptly described by a "cost function." We define probability 
distributions of the two directional parameters, called generating distributions since they generate possible 
valleys or states we are to explore. We define another distribution, called the acceptance distribution, 
which depends on the difference of cost functions of the present generated valley we are to explore and 
the last saved lowest valley. The acceptance distribution decides probabilistically whether to stay in a new 
lower valley or to bounce out of it. All the generating and acceptance distributions depend on 
temperatures. 

In 1984 [47], it was established that SA possessed a proof that, by carefully controlling the rates of 
cooling of temperatures, it could statistically find the best minimum, e.g., the lowest valley of our 
example above. This was good news for people trying to solve hard problems which could not be solved 
by other algorithms. The bad news was that the guarantee was only good if they were willing to run SA 
forever. In 1987, a method of fast annealing (FA) was developed [48], which permitted lowering the 
temperature exponentially faster, thereby statistically guaranteeing that the minimum could be found in 
some finite time. However, that time still could be quite long. Shortly thereafter, Very Fast Simulated 
Reannealing (VFSR) was developed in 1987 [45], now called Adaptive Simulated Annealing (ASA), 
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which is exponentially faster than FA. 

ASA has been applied to many problems by many people in many disciplines [49-51]. The 
feedback of many users regularly scrutinizing the source code ensures its soundness as it becomes more 
flexible and powerful. 

5.1.2. Mathematical outline 

ASA considers a parameter a\ in dimension i generated at annealing-time k with the range 

aie[Ai,Bi], (30) 
calculated with the random variable y l , 

y'e[-l,l]. (31) 

The generating function griy) is defined, 

d i D 

(=1 2(ly ; l + T i )ln(l + l/r ! ) M 

where the subscript i on T t specifies the parameter index, and the ^-dependence in Tj(k) for the annealing 
schedule has been dropped for brevity. Its cumulative probability distribution is 



G T (y) = J • • • J dy' 1 ■ ■ ■ dy' D g T (y) =n<W) , 



-l -l 

Gk y ) = l + ^'>^™> . (33) 

2 2 ln(l + l/r ; ) 
y' is generated from a u l from the uniform distribution 

u'eU[0,l] , 

y = sgn (V - i)T,.[(l + l/r,) 12 "'" 1 ' - 1] • (34) 

It is straightforward to calculate that for an annealing schedule for T t 
T i (k) = T 0ie xp(-c i k l/D ), (35) 
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a global minima statistically can be obtained. I.e., 



oo oo 



z^-Z[n^V ] ^ =o °- (36) 

k k Q i=\ 2\y l \ci k 
Control can be taken over c u such that 

T fi = T Qi exp(-m ; ) when k f = exp n ; , 

Ci = mi exp(-« ( /D) , (37) 
where m ; and n t can be considered "free" parameters to help tune ASA for specific problems. 

ASA has over 100 OPTIONS available for tuning. A few important ones were used in this project. 

5.1.3. Reannealing 

Whenever doing a multi-dimensional search in the course of a complex nonlinear physical problem, 
inevitably one must deal with different changing sensitivities of the a' in the search. At any given 
annealing-time, the range over which the relatively insensitive parameters are being searched can be 
"stretched out" relative to the ranges of the more sensitive parameters. This can be accomplished by 
periodically rescaling the annealing-time k, essentially reannealing, every hundred or so acceptance- 
events (or at some user-defined modulus of the number of accepted or generated states), in terms of the 
sensitivities s t calculated at the most current minimum value of the cost function, C, 

Si = dC/da* . (38) 

In terms of the largest s t = s max , a default rescaling is performed for each k { of each parameter dimension, 
whereby a new index k' t is calculated from each k t , 

kj > k j , 

T ik' — T ifc^S max l S i) , 

k'i = (ln(T i{) /T ik ,)/ Ci ) D . (39) 
r,o is set to unity to begin the search, which is ample to span each parameter dimension. 
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5.1.4. Quenching 

Another adaptive feature of ASA is its ability to perform quenching in a methodical fashion. This 
is applied by noting that the temperature schedule above can be redefined as 

T i (k i ) = T 0ie xp(-c i kf-' D ), 

c i = m i exp(-n i Q i /D) , (40) 
in terms of the "quenching factor" Q t . The sampling proof fails if <2; > 1 as 

z n i/£ e,/D = e i/£ a < oo . (4i) 

k k 

This simple calculation shows how the "curse of dimensionality" arises, and also gives a possible 
way of living with this disease. In ASA, the influence of large dimensions becomes clearly focussed on 
the exponential of the power of k being 1/D, as the annealing required to properly sample the space 
becomes prohibitively slow. So, if resources cannot be committed to properly sample the space, then for 
some systems perhaps the next best procedure may be to turn on quenching, whereby Q t can become on 
the order of the size of number of dimensions. 

The scale of the power of 1/D temperature schedule used for the acceptance function can be altered 
in a similar fashion. However, this does not affect the annealing proof of ASA, and so this may used 
without damaging the sampling property. 

5.2. x-Indicator of Market Contexts 

Our studies of contexts of markets well recognized by option traders to have significantly different 
volatility behavior show that the exponents x are reasonably faithful indicators defining these different 
contexts. 

We feel the two-factor model is more accurate because the data indeed demonstrate stochastic 
volatility [4]. We also note that the two-factor x's are quite robust and uniform when being fit by ASA 
across the last few years. This is not true of the one-factor ASA fitted x's unless we do not use the Black- 
Scholes a as a parameter, but rather calculate as historical volatility during all runs. Some results of two- 
factor studies and one-factor studies using a Black-Scholes a have been reported elsewhere [13]. 
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Since a is not widely traded and arbitraged, to fit the two-factor model, we calculate this quantity 
as an historical volatility for both its prepoint and postpoint values. Some previous studies used a scaled 
implied volatility (which is calculated from a BS model). We use a standard deviation o' , 

a' = StdDev(dS / F(S, S , S^, x, y)) . (42) 

In the one-factor model, it does not make good numerical sense to have two free parameters in one 
term, i.e., a and x, as these cannot be fit very well within the variance the data. Instead, one method is to 
take guidance from the two-factor results, to set a scale for an effective a, and then fit the parameter x. 
Another method it apply the above StdDev as a proxy for a. Some motivation for this approach is given 
by considering collapsing a two-factor stochastic volatility model in one-factor model: The one-factor 
model now has an integral over the stochastic process in its diffusion term. The is integral is what we are 
approximating by using the standard deviation of a moving window of the data. 

6. PATH-INTEGRAL (PATHINT) DEVELOPMENT 
6.1. PATHINT Outline 

The fits described above clearly demonstrate the need to incorporate stochastic volatility in option 
pricing models. If one-factor fits are desired, e.g., for efficiency of calculation, then at the least the 
exponent of price x should be permitted to freely adapt to the data. In either case, it is required to develop 
a full set of Greeks for trading. To meet these needs, we have used a path-integral code, PATHINT, 
described below, with great success. At this time, the two-factor code takes too long to run for daily use, 
but it proves to be a good weekly baseline for the one-factor code. 

The PATHINT algorithm develops the long-time probability distribution from the Lagrangian fit by 
the first optimization code. A robust and accurate histogram-based (non-Monte Carlo) path-integral 
algorithm to calculate the long-time probability distribution has been developed to handle nonlinear 
Lagrangians [18-20,22,52-54], 

The histogram procedure recognizes that the distribution can be numerically approximated to a high 
degree of accuracy as sum of rectangles at points M i of height P, and width AM,-. For convenience, just 
consider a one-dimensional system. The above path-integral representation can be rewritten, for each of 
its intermediate integrals, as 
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P(M; t + AO = J dM> [g l-(2xAtT m exp(-L^)] W; t) = J dM>G(M, M'\ At)P(M>; t) , 



P(M;t) = jt*W-M i )P i (t), 



n(M - Md = ■ 



, (M; - - AM,, ) < M < (M,- + - AM,) , 

2 2 , (43) 

1 , otherwise , 



which yields 

P i (t + At) = T ij (At)P j (t), 

TJAt) = f M,+AM,/2 dM f Mj+AMj/2 dM ' G ( M ^ M >- Af ) _ ( 44 ) 

1 AM;.] + AM; J M—AMj_i/2 J Mj-AM^Il 

Tjj is a banded matrix representing the Gaussian nature of the short-time probability centered about the 
(varying) drift. 

Care must be used in developing the mesh in AM G , which is strongly dependent on the diagonal 
elements of the diffusion matrix, e.g., 

AM G « (A^ IGIIGI ) 1/2 . (45) 

Presently, this constrains the dependence of the covariance of each variable to be a nonlinear function of 
that variable, albeit arbitrarily nonlinear, in order to present a straightforward rectangular underlying 
mesh. Below we address how we have handled this problem in our two-factor stochastic-volatility model. 

Fitting data with the short-time probability distribution, effectively using an integral over this 
epoch, permits the use of coarser meshes than the corresponding stochastic differential equation. The 
coarser resolution is appropriate, typically required, for numerical solution of the time-dependent path- 
integral: By considering the contributions to the first and second moments of AM G for small time slices 9, 
conditions on the time and variable meshes can be derived [52]. The time slice essentially is determined 
by 6 < L~ l , where L is the "static" Lagrangian with dM G ldt = 0, throughout the ranges of M G giving the 
most important contributions to the probability distribution P. The variable mesh, a function of M G , is 
optimally chosen such that AM is measured by the covariance g , or AM ~(g 0) . 
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If the histogram method is further developed into a trapezoidal integration, then more accuracy can 
be expected at more complex boundaries [53]. Such problems does not arise here, and 6-7 significant 
figure accuracy is easily achieved provided great care is taken to develop the mesh as described above. 
For example, after setting the initial-condition discretized delta function MSUPG at the prepoint of an 
interval, the mesh going forward in M G is simply calculated stepwise using 

AM G = (g GG 0) m . (46) 

However, going backwards in M G , an iterative procedure was used at each step, starting with an estimate 
from the prepoint and going forward again, until there was no mismatch. That much care is required for 
the mesh was observed in the original Wehner-Wolfer paper [52]. 

It is important to stress that very good numerical accuracy is required to get very good Greeks 
required for real-world trading. Many authors develop very efficient numerical schemes to get reasonable 
prices to 2 or 3 significant figures, but these methods often are not very good to enough significant figures 
to get good precision for the Greeks. Typical Monte Carlo methods are notorious for giving such poor 
results after very long computer runs. In particular, we do not believe that good Greeks required for 
trading can be obtained by using meshes obtained by other simpler algorithms [55]. 

The PATHINT algorithm in its present form can "theoretically" handle any n-factor model subject 
to its diffusion-mesh constraints. In practice, the calculation of 3-factor and 4-factor models likely will 
wait until giga-hertz speeds and giga-byte RAM are commonplace. 

6.2. Development of Long-Time Probabilities 

The noise determined empirically as the diffusion of the data is the same, independent of x within 
our approach. Therefore, we scale different exponents such that the the diffusions, the square of the 
"basis-point volatilities" (BPV), are scaled to be equivalent. Then, there is not a very drastic change in 
option prices for different exponents x for the strike X set to the 5 underlying, the at-the-money (ATM) 
strike. This is not the case for out of the money (OTM) or in the money (ITM) strikes, e.g., when 
exercising the strike would generate loss or profit, resp. This implies that current pricing models are not 
radically mispricing the markets, but there still are significant changes in Greeks using more sophisticated 
models. 
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6.3. Dependence of Probabilities on 5 and x 

Fig. 2 gives examples of the short-time distribution evolved out to T = 0. 5 year for x in {-1, 0, 1, 
2}, with 500 intermediate epochs/foldings, and BS a = 0. 0075. Each calculation scales a by multiplying 
by S/F(S, So, S^, x, y). 



Fig. 2. 



Fig. 3 gives an example of a two-factor distribution evolved out to T = 0. 5 year for x = 0. 7. 



Fig. 3. 



6.4. Two-Factor Volatility and PATHINT Modifications 

In our two-factor model, the mesh of 5 would depend on o and cause some problems in any 
PATHINT grid to be developed in S-a. 

For some time we have considered how to handle this generic problem for rc-factor multivariate 
systems with truly multivariate diffusions within the framework of PATHINT. In one case, we have taken 
advantage of the Riemannian invariance of the probability distribution as discussed above, to transform to 
a system where the diffusions have only "diagonal" multiplicative dependence [19]. However, this leads 
to cumbersome numerical problems with the transformed boundary conditions [20]. Another method, not 
yet fully tested, is to develop a tiling of diagonal meshes for each factor i that often are suitable for off- 
diagonal regions in an w-factor system, e.g., 

AM[ = 2 m *AM , 

AM » ^gf»At , (47) 

where the mesh of variable i at a given point labeled by k is an exponentiation of 2, labeled by m\\ the 
integral power m\ is determined such that it gives a good approximation to the diagonal mesh given by 
the one-factor PATHINT mesh conditions, in terms of some minimal mesh AM' , throughout regions of 
the Lagrangian giving most important contributions to the distribution as predetermined by a scan of the 
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system. This tiling of the kernel is to be used together with interpolation of intermediate distributions. 

The results of our study here are that, after the at-the-money BPV are scaled to be equivalent, there 
is not a very drastic change in the one-factor ATM Greeks developed below. Therefore, while we have 
not at all changed the functional dependence of the Lagrangian on S and a, we have determined our 
meshes using a diffusion for the 5 equation as cr F(S, Sq, S^, x, y), where a is determined by the same 
BPV-equivalent condition as imposed on the one-factor models. This seems to work very well, especially 
since we have taken our cr equation to be normal with a limited range of influence in the calculations. 
Future work yet has to establish a more definitive distribution for cr. 

7. CALCULATION OF DERIVATIVES 

7.1. Primary Use of Probabilities For European Options 

We can use PATHINT to develop the distribution of the option value back in time from expiration. 
This is the standard approach used by CRR, explicit and implicit Crank-Nicolson models, etc [56]. 

For European options, we also take advantage of the accuracy of PATHINT enhanced by 
normalizing the distribution as well as the kernel at each iteration (though in these studies this was not 
required after normalizing the kernel). Therefore, we have calculated our European option prices and 
Greeks using the most elementary and intuitive definition of the option's price V [57], which is the 
expected value 



where X is the strike price, and the expected value < . > is taken with respect to the risk-neutral 
distribution of the underlying market 5. It should be noted that, while the standard approach of 
developing the option price delivers at the present time a range of underlying values for a given strike, our 
approach delivers a more practical complete range of strikes (as many as 50-60 for Eurodollar options) 
for a given underlying at the present time, resulting in a greatly enhanced numerical efficiency. The risk- 
neutral distribution is effectively calculated taking the drift as the cost-of-carry b times S, using the above 
arguments leading to the BS formula. We have designed our codes to use parameters risk-free -rate r and 
cost-of-carry b such that 




put 



call 



(48) 
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b = r, cost of carry on nondividend stock , 

b = r-q, cost of carry on stock paying dividend yield q , 

b = 0, cost of carry on future contract , 

b = r-rf, cost of carry on currency with rate r y , (49) 
which is similar to how generalized European BS codes use b and r [58]. 
Using this approach, the European price V E is calculated as 

V E =< max[z(e (b - r)T S T - e~ rT X), 0] > . (50) 

The American price V A must be calculated using a different kernel going back in time from 
expiration, using as "initial conditions" the option values used in the above average. This kernel is the 
transposed matrix used for the European case, and includes additional drift and "potential" terms due to 
the need to develop this back in time. This can be understood as requiring the adjoint partial differential 
equation or a postpoint Lagrangian in real time. 

That is, a forward equation for the conditional probability distribution P[M(t + dt), t + dt\M{t), t] is 

dP/dt = l - (g GG >), GG < - (g G P\ G + V, (51) 

where the partial derivatives with respect to M G act on the postpoint M G (t + dt). A backward equation 
for the conditional probability distribution P[M(t + dt), t + dt\M(t), t] is 

dP/dt = l - g GG > GG - - g G P G + V , (52) 

where the partial derivatives with respect to M act on the prepoint M (t). The forward equation has a 
particularly simple form in the mathematically equivalent prepoint path integral representation. 

Above, we have described how the forward distribution at present time T Q is evolved using 
PATHINT to the time of expiration, P(T), e.g., using a path-integral kernel K folded over n epochs, 
where it is folded with the function O, 

V=< 0(T) P(T) > = < O(K n P(t )) > , 

0(T) = max[z(e (b - r)T S T - e~ rT X), 0] , (53) 
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to determine the European value at the present time of the calls and puts at different strike values X. An 
equivalent calculation can be performed by using the backward equation, expressed in terms of the 
"equivalent" kernel acting on O, 



It is convenient to use the simple prepoint representation for the Lagrangian, so the backward kernel is 
first re-expressed as a forward kernel by bringing the diffusions and drifts inside the partial derivatives, 
giving a transformed adjoint kernel K\ 

The above mathematics is easily tested by calculating European options going forwards and 
backwards. For American options, while performing the backwards calculation, at each point of mesh, 
the options price evolving backwards from T is tested and calculated as 



The Greeks {A, F, 0} are directly taken off the final developed option at the present time, since the M G 
mesh is available for all derivatives. We get excellent results for all Greeks. Note that for CRR trees, 
only one point of mesh is at the present time, so A requires moving back one epoch and F requires moving 
back two epochs, unless the present time is pushed back additional epochs, etc. 

7.2. PATHINT Baselined to CRR and BS 

The CRR method is a simple binomial tree which in a specific limit approaches the BS partial 
differential equation. It has the virtues of being fast and readily accommodates European and American 
calculations. However, it suffers a number of well-known numerical problems, e.g., a systematic bias in 
the tree approximation and an oscillatory error as a function of the number of intermediate 
epochs/iterations in its time mesh. Some Greeks like {A, F, 0} can be directly taken off the tree used for 
pricing with reasonable approximations (at epochs just before the actual current time). The first problem 
for American options can be alleviated somewhat by using the variant method [5], 

CRRyariant = CRRAmerican _ CRR European + BS . (56) 

The second problem can be alleviated somewhat by averaging runs of n iterations with runs of n + 1 
iterations [59]. This four-fold increase of runs is rarely used, though perhaps it should be more often. 
Furthermore, if increased accuracy in price is needed in order to take numerical derivatives, typically 



V =< 0(T) P(T) > = < (K n O(T))P(t ) > . 



(54) 



O new = max[(S-X),(e- Mf O old )] . 



(55) 
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200-300 iterations should be used for expirations some months away, not 30-70 as too often used in 
practice. 

When taking numerical derivatives there can arise a need to tune increments taken for the 
differentials. For some Greeks like A and T the size of the best differentials to use may vary with strikes 
that have different sensitivities to parts of the underlying distribution. One method of building in some 
adaptive flexibility across many such strikes is to increase the order of terms used in taking numerical 
derivatives. (This was not required for results presented here.) For example, it is straightforward to verify 
that, while the central difference 

df_ = f{x + dx) - fjx - dx) 
dx 2dx 

is typically good to o((<ix) 3 ), 

df _ -fix + 2dx) + 8 fix + dx) - 8 fix - dx) + fix - 2dx) 
dx \2dx 

is typically good to o(idx) 5 ). Similarly, while 

d 2 f = fjx + dx)-2fix) + fjx-dx) 
dx 2 dx dx 

is typically good to ((c?x) 4 ), 

d 2 f -dix + 2dx) + 16 fix + dx) - 30 fix) + 16 fix - dx) - fix - 2dx) 

(60) 



dx 2 dx dx 



is typically good to Hdx) ). 



Table 1 gives an example of baselining our one-factor PATHINT code to the CRR and BS codes 
using the above safeguards for an option whose American price is the same as its European counterpart, a 
typical circumstance [58]. In the literature, the CRR code is most often taken as the benchmark for 
American calculations. We took the number of intermediate epochs/points to be 300 for each calculation. 
Parameters used for this particular ATM call are T = 0. 5 years, r = 0. 05, b = 0, a = 0. 10. 



Table 1. 
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Tests with American CRR and American PATHINT lead to results with the same degrees of 
accuracy. 

7.3. Two-Factor PATHINT Baselined to One-Factor PATHINT 

Previous papers and tests have demonstrated that the two-factor PATHINT code performs as 
expected. The code was developed with only a few lines to be changed for running any ^-factor problem, 
i.e., of course after coding the Lagrangian specific to a given system. Tests were performed by combining 
two one-factor problems, and there is no loss of accuracy. However, here we are making some additional 
mesh approximations as discussed above to accommodate a in the 5 diffusion. This seems quite 
reasonable, but there is no sure test of the accuracy. We indeed see that the ATM results are very close 
across x's, similar to our ATM comparisons between BS and our one-factor PATHINT results for various 
x's, where again scaling is performed to have all models used the same BPV (using the cr procedure for 
the mesh as described above for the two-factor model). 

The logical extension of Greeks for the two-factor model is to develop derivatives of price with 
respect to p and e in a volatility equation. However, we did not find a bona fide two-factor proxy for the 
one-factor T, the derivative of price with respect to the one-factor a constant. We get very good ATM T 
comparisons between BS and our one-factor models with various x's. We tried simply multiplying the 
noise in the two-factor stochastic volatility in the price equation by a parameter with deviations from 1 to 
get numerical derivatives of PATHINT solutions, and this gave somewhat good agreement with the ATM 
BPV-scaled BS T within a couple of significant figures. Perhaps this is not too surprising, especially 
given the correlation substantial p between the price and volatility equations which we do not neglect. 

8. CONCLUSIONS 

The results of our study are that, after the at-the-money basis-point volatilities are scaled to be 
equivalent, there is only a very small change in option prices for different exponents x, both for the one- 
factor and two-factor models. There still are significant differences in Greeks using more sophisticated 
models, especially for out-of-the-money options. This implies that current pricing models are not 
radically mispricing the markets. 
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Our studies point to contexts of markets well recognized by option traders to have significantly 
different volatility behavior. Suppression of stochastic volatility in the one-factor model just leaks out 
into stochasticity of parameters in the model, e.g., especially in x, unless additional numerical methods 
are employed, e.g., using an adaptive standard deviation. Our studies show that the two-factor exponents 
x are reasonably faithful indicators defining these different contexts. The x -exponents in the two-factor 
fits are quite stable. As such, especially the two-factor x can be considered as a "context indicator" over a 
longer time scale than other indicators typically used by traders. The two-factor fits also exhibit 
differences due to the a parameters, including the p correlations, in accord with the sense traders have 
about the nature of changing volatilities across this time period. 

Modern methods of developing multivariate nonlinear multiplicative Gaussian-Markovian systems 
are quite important, as there are many such systems and the mathematics must be diligently exercised if 
such models are to faithfully represent the data they describe. Similarly, sophisticated numerical 
techniques, e.g., global optimization and path integration are important tools to carry out the modeling 
and fitting to data without compromising the model, e.g., by unwarranted quasi-linear approximations. 
Three quite different systems have benefited from this approach: 

The large-scale modeling of neocortical interactions has benefited from the use of intuitive 
constructs that yet are faithful to the complex algebra describing this multiple- scaled complex system. 
For example, canonical-momenta indicators have been successfully applied to multivariate financial 
markets. 

It is clear that ASA optimization and PATHINT path-integral tools are very useful to develop the 
algebra of statistical mechanics for a large class of nonlinear stochastic systems encountered in finance. 
However, it also is clear that each system typically presents its own non-typical unique character and this 
must be included in any such analysis. A virtue of this statistical mechanics approach and these 
associated tools is they appear to be flexible and robust to handle quite different systems. 
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FIGURE CAPTIONS 

FIG. 1. (a) F(S, Sq, S^, x > y) for x = 1, the Black-Scholes case. The other parameters are 5 = 5, 
S = 0. 5, 5^ = 20, y = 0. (b) F(S, S , S^, x, y) for x = 0, the normal distribution, (c) F(S, S , S^, x, y) 
for x = —l. (d) F(S, Sq, S^, x, y) for x = 2. 

FIG. 2. The short-time probability distribution at time T = 0.5 years for x= 1, the (truncated) 
Black-Scholes distribution. The short-time probability distribution at time T = 0. 5 years for x = 0, the 
normal distribution. The short-time probability distribution at time T = 0. 5 years for x = -1. The short- 
time probability distribution at time T = 0. 5 years for x = 2. 

FIG. 3. A two-factor distribution evolved out to T = 0. 5 year for x = 0. 7. 
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TABLE CAPTIONS 

Table 1. Calculation of prices and Greeks are given for closed form BS (only valid for European 
options), binomial tree CRR European , CRR Amer i can , CRR var ; ant , and PATHINT. As verified by calculation, 
the American option would not be exercised early, so the PATHINT results are identical to the European 
option. The CRR Amer i C an differs somewhat from the CRR European due to the discrete nature of the 
calculation. All CRR calculations include averaging over 300 and 301 iterations to minimize oscillatory 
errors. 



High-resolution path-integral . 



- Figure 1 - 



Lester Ingber 




High-resolution path-integral . 



- Figure 2 - 



Lester Ingber 




High-resolution path-integral . 



- Figure 3 - 



Lester Ingber 



Two-Factor Probability 



Long-Time Probability 



100 r 
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Greek 


BS 


CRPvEuropean 


CRRAmerican 


CRRvariant 


PATHINT 


Price 


0.138 


0.138 


0.138 


0.138 


0.138 


Delta 


0.501 


0.530 


0.534 


0.506 


0.501 


VJalllllla 


1 100 


1 14? 


1 1 SQ 


1 1 1fS 


1 100 


Theta 


-0.131 


-0.130 


-0.132 


-0.133 


-0.131 


Rho 


-0.0688 


-0.0688 


-0.0530 


-0.0530 


-0.0688 


Vega 


1.375 


1.375 


1.382 


1.382 


1.375 



